Polarization effects in protein-ligand calculations extend farther than the actual induction energy

The various roles that polarizabilities play in the calculation of protein-ligand interaction energies with a polarizable force field are investigated, and the importance and distance dependence of some common approximations is determined for each of these roles separately, using quantum-mechanical calculations as the reference. It is found that the pure induction energy, if defined as the energetic gain from the charge redistribution upon interaction between the protein and ligand, is a rather short-ranged effect that becomes independent of the exact implementation at distances above ∼4Å. On the other hand, the polarization between the protein residues in the assembly of a protein from separately computed fragments (as is routinely done in force field development) has a significant effect on the computed interaction energies, even for residues as far as 15Å from the ligand. Finally, polarization improves the transferability of partial charges, but the simple polarization model used in, for example, the Amber force field explains only 14-19% of the conformational variation of the charges. In all cases, more advanced polarization models, especially involving anisotropic polarizabilities, seem to give significantly better descriptions of these effects. The study suggests that an accurate treatment of polarization can be important even in systems where the actual induction energy is small in magnitud

Accurate reaction energies in proteins obtained by combining QM/MM and large QM calculations

We here suggest and test a new method to obtain stable energies in proteins for charge-neutral reactions by running large quantum mechanical (QM) calculations on structures obtained by combined QM and molecular mechanics (QM/MM) geometry optimisation on several snapshots from molecular dynamics simulations. As a test case, we use a proton transfer between a metal-bound cysteine residue and a second-sphere histidine residue in the active site of [Ni,Fe] hydrogenase, which has been shown to be very sensitive to the surroundings. We include in the QM calculations all residues within 4.5 Å of the active site, two capped residues on each side of the active-site residues, as well as all charged groups that are buried inside the protein, which for this enzyme includes three iron–sulphur clusters, in total 930 atoms. These calculations are performed at the BP86/def2-SV(P) level, but the energies are then extrapolated to the B3LYP/def2-TZVP level with a smaller QM system and zero-point energy, entropy, and thermal effects are added. We test three approaches to model the remaining atoms of the protein solvent, viz. by standard QM/MM approaches using either mechanical or electrostatic embedding, or by using a continuum solvation model for the large QM systems. Quite encouragingly, the three approaches give the same results within 13 kJ/mol and variations in the size of the QM system do not change the energies by more than 8 kJ/mol, provided that the QM/MM junctions are not moved closer to the QM system. The statistical precision for the average over ten snapshots is 1–3 kJ/mol

Ligand affinities estimated by quantum chemical calculations

We present quantum chemical estimates of ligand-binding affinities performed, for the first time, at a level of theory for which there is a hope that dispersion and polarization effects are properly accounted for (MP2/cc-pVTZ) and at the same time effects of solvation, entropy, and sampling are included. We have studied the binding of seven biotin analogues to the avidin tetramer. The calculations have been performed by the recently developed PMISP approach (polarizable multipole interactions with supermolecular pairs), which treats electrostatic interactions by multipoles up to quadrupoles, induction by anisotropic polarizabilities, and nonclassical interactions (dispersion, exchange repulsion, etc.) by explicit quantum chemical calculations, using a fragmentation approach, except for long-range interactions that are treated by standard molecular-mechanics Lennard-Jones terms. In order to include effects of sampling, 10 snapshots from a molecular dynamics simulation are studied for each biotin analogue. Solvation energies are estimated by the polarized continuum model (PCM), coupled to the multipole-polarizability model. Entropy effects are estimated from vibrational frequencies, calculated at the molecular mechanics level. We encounter several problems, not previously discussed, illustrating that we are first to apply such a method. For example, the PCM model is, in the present implementation, questionable for large molecules, owing to the use of a surface definition that gives numerous small cavities in a protein

Effect of geometry optimisations on QM-cluster and QM/MM studies of reaction energies in proteins

We have examined the effect of geometry optimisation on energies calculated with the quantum-mechanical (QM) cluster, the combined QM and molecular-mechanics (QM/MM), the big-QM approaches (very large single-point QM calculations taken from QM/MM-optimised structures, including all atoms within 4.5 Å of the minimal active site, all buried charged groups in the protein, and truncations moved at least three residues away from the active site). We study a simple proton-transfer reaction between His-79 and Cys-546 in the active site of [Ni,Fe] hydrogenase and optimise QM systems of 50 different sizes (56–362 atoms). Geometries optimised with QM/MM are stable and reliable, whereas QM-cluster optimisations give larger changes in the structures and sometimes lead to large distortions in the active site if some hydrogen-bond partners to the metal ligands are omitted. Keeping 2–3 atoms for each truncated residue (rather than one) fixed in the optimisation improves the results, but does not solve all problems for the QM-cluster optimisations. QM-cluster energies in vacuum and a continuum solvent are insensitive to the geometry optimisations with a mean absolute change upon the optimisations of only 4–7 kJ/mol. This shows that geometry optimisations do not decrease the dependence of QM-cluster energies on how the QM system is selected – there is still a ~60 kJ/mol difference between calculations in which groups have been added to the QM system according to their distance to the active site or based on QM/MM free-energy components. QM/MM energies do not show such a difference, but they converge rather slowly with respect to the size of the QM system, although the convergence is improved by moving truncations away from the active site. The big-QM energies are stable over the 50 different optimised structures, 57±1 kJ/mol, although some smaller trends can be discerned. This shows that both QM-cluster geometries and energies should be interpreted with caution. Instead, we recommend QM/MM for geometry optimisations and energies calculated by the big-QM approach

Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model.

The most general way to improve the accuracy of binding-affinity calculations for protein-ligand systems is to use quantum-mechanical (QM) methods together with rigorous alchemical-perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin-3 at a reasonably high QM level (dispersion-corrected density functional theory with a triple-zeta basis set) and with a sufficiently large QM system to include all short-range interactions with the ligand (744-748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non-Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ∼3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations. © 2015 Wiley Periodicals, Inc

Conformational Dependence of Isotropic Polarizabilities

We perform a statistical and energetic analysis of atomic polarizabilities obtained with the LoProp approach for all atoms in the avidin tetramer for 70 snapshots from molecular dynamics simulations with seven different biotin analogues, and from the crystal structure of the photosynthetic reaction center (in total 560 698 individual polarizabilities). Dynamic effects give a variation of the polarizabilities of 0.09 angstrom(3) on average. Atoms at different positions in the sequence show a variation of 0.14 angstrom(3) on average, caused by the conformational dependence of the polarizabilities. This variation gives errors of 2 and 1 kJ/mol for relative conformational and ligand-binding induction energies. Averaged elementwise or atom type polarizabilities give larger errors, e.g., 9 and 7 kJ/mol, respectively, for the relative conformational energies. Therefore, we recommend that polarizabilities should be assigned atomwise (i.e., individual polarizabilities for each atom in all residues), in the same way as for charges. We provide such a set of extensively averaged polarizabilities (xAvPol) for all atoms in avidin and the photosynthetic reaction center, applicable at the B3LYP/aug-cc-pVTZ level, which is converged with respect to the basis-set limit

Converging ligand-binding free energies obtained with free-energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free-energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa-acid deep-cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158-224 atoms). We use single-step exponential averaging (ssEA) and the non-Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi-empirical PM6-DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free-energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions

A first-principles approach to protein–ligand interaction

It is still impossible to make an accurate, purely theoretical prediction of the free energy of a ligand binding to a protein in aqueous environment. The two main problems are the immense number of nuclear configurations contributing to the binding free energy and the impossibility to apply accurate quantum-chemical methods to such a large system, even for a single configuration. In this thesis, the second of these problems is addressed by exploring various ways of approximating the quantum-chemical interaction energy without introducing experimental data in the models. The use of quantum chemistry to derive parameters for advanced molecular mechanics models is explored. First, models for the repulsion term, based on either orbital overlap or electron density overlap, are compared. The latter models are found to be inaccurate for certain interactions, although they perform well on average. Second, the distributed multipoles and polarizabilities required for the electrostatic and induction terms are assessed, the result suggesting that the newly developed LoProp method is an improvement on earlier methods. The accuracy of various approximations inherent in the polarizability model are also tested. It is found that the neglect of Pauli effects is a severe approximation, but that the polarizability model nevertheless gives reasonable results, owing to error cancellation. As a basis for future polarization models that avoid such cancellation, a quantum-chemical model is introduced, in which Pauli effects from surrounding molecules are included through a pseudopotential. Finally, a method for protein–ligand interactions is developed, in which the protein is divided into fragments and the pair potentials between the ligand and each fragment is calculated by quantum chemistry, whereas the non-additivity is modeled by multipoles and polarizabilities. This and further approximations are tested, allowing for a full protein–ligand interaction energy to be computed at an unprecedented level of theory. The method is applied in an approximate calculation of binding free energies for a set of ligands to avidin, unfortunately giving poor results. A possible reason for this failure is the treatment of solvation